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Abstract 

Max-stable processes have proved to be useful for the statistical modelling of spatial 
extremes. Several representations of max-stable random fields have been proposed in the 
literature. One such representation is based on a limit of normalized and scaled pointwise 
maxima of stationary Gaussian processes that was first introduced by Kabluchko, Schlather 
and de Haan [Mi . 

This paper deals with statistical inference for max-stable space-time processes that are de- 
fined in an analogous fashion. We describe pairwise likelihood estimation, where the pairwise 
density of the process is used to estimate the model parameters and prove strong consistency 
and asymptotic normality of the parameter estimates for an increasing space-time dimension, 
i.e., as the joint number of spatial locations and time points tends to infinity. A simulation 
study shows that the proposed method works well for these models. 
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1 Introduction 



Max-stable processes have proven to be useful in the modelling of spatial extremes. Typically, 
meteorological extremes like heavy rainfall or extreme wind speeds are modelled using extreme 
value theory. In particular, sample maxima such as annual maximum wind speeds are observed 
at several locations of some spatial process. Other applications may involve the analysis of image 
data resulting from tomographic examinations. 

Several representations of max-stable processes have been proposed in the literature, includ- 
ing for example Brown and Resnick [3j, de Haan [SJ, Kabluchko, Schlather and de Haan 
and Schlather |19| . Recently, models for extreme values observed in a space-time setting have 
generated a great deal of interest. First approaches can be found in Davis and Mikosch [6], Huser 
and Davison jllj . Kabluchko [T3], and Davis, Kliippelberg and Steinkohl [5]. 

In this paper, we follow the approach described in Davis et al. [5], who extend the max-stable 
process introduced in Kabluchko et al. |14) to a space-time setting. The process is constructed 
as the limit of rescaled and normalized maxima of independent replications of some stationary 
Gaussian space-time process. The underlying correlation function of the Gaussian process is as- 
sumed to belong to a parametric model whose parameters describe smoothness of the correlation 
function near the origin. 

As it is well-known for max-stable processes, the full likelihood function is computationally 
intractable and other methods have to be used to derive parameter estimates. Standard pro- 
cedures for such cases are composite likelihood including pairwise likelihood estimation. These 
methods go back to Besag [1], and there is an extensive literature available dealing with applica- 
tions and properties of the estimates, see for example Gox and Reid [1], Lindsay [16], Varin [23] . 
or Varin and Vidoni [23] . Recent work concerning the application of pairwise likelihood methods 
to max-stable random fields can be found in Huser and Davison ^llj and Padoan, Ribatet and 
Sisson [18j. 

Since the observations in a space-time setting are correlated in space and time, we use special 
properties of max-stable processes to show strong consistency and asymptotic normality of the 
estimates. Here, it is assumed that the locations lie on a regular lattice and that the time points 
are equidistant. The spatial and/or the temporal dimension, i.e., the number of spatial locations 
and/or time points, increases to infinity. The main step of the proof is based on a strong law 
of large numbers for the pairwise likelihood function. Stoev |20j analyzed ergodic properties 
for max-stable processes in time resulting from extremal integral representations for max-stable 
processes that were introduced in Stoev and Taqqu [21] . The extension to a spatial setting and 
the resulting strong law of large numbers was shown by Wang, Roy and Stoev [26]. By combining 
these two results we obtain a strong law of large numbers for a jointly increasing space-time 
domain. 

In addition to strong consistency, we prove asymptotic normality for the pairwise likelihood 
estimates. A first result concerning asymptotic normality of pairwise likelihood estimates for 
max-stable space-time processes can be found Huser and Davison fTT], who fix the number of 
locations and let the number of time points tend to infinity. We formulate asymptotic normality 
for the space-time setting and use Bolthausen's theorem [2j together with strong mixing prop- 
erties shown by Dombry and Eyi-Minko [9] to prove asymptotic normality for an increasing 
number of space-time locations. 

Our paper is organized as follows. In Section [2| we introduce the max-stable space-time 
process for which inference properties will be considered in subsequent sections. Section |3.2| 
describes pairwise likelihood estimation and the particular setting for our model. In Section [4] 
we establish strong consistency for the estimates for increasing space-time domain. Asymptotic 
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normality of these parameters is established in Section [5| A simulation study evaluating the 
performance of the estimates is described in Section [6| 



2 Description of the model 

We start with the process that will be used for modelling extremes in space and time; details 
can be found in Davis et al. [5]. Let {Z{s,t),s G M.'^,t G [0,oo)} denote a stationary space-time 
Gaussian process on x [0, oo) with mean zero and variance one. With the correlation function 

p{h, u)=E [Z{s, t)Z{s + h,t + u)], 

where /i G M'^ is the spatial lag and u G M is the time lag, we make the following assumption 
that will be used throughout the paper. 

Assumption 2.1. There exist sequences of constants s„ — )• 0, t„ — )• as n ^ oo, such that 

log n(l — p{snh, tnu)) — )• 5{h, u) > 0, as n ^ oo. 



Assumption 2.1 is natural in the context of stationary space-time models; the correlation 
function tends to one at a certain rate as the space-time lag approaches the zero. 

Proposition 2.2 (Kabluchko et al. pjj and Davis et al. [5]). Let {Zj{s, t),s eW^,te [0, oo)} , j = 
1, . . . ,n, he independent replications of the space-time Gaussian process described above and let 
£ N} denote points of a Poisson random measure on [0, oo) with intensity measure ^'"^d^. 
Suppose Assumption 2.1 is satisfied. Then, the random fields {r]n{s,t)), s £R'^,tG [0,oo)}, de- 
fined for n £N by 

n ^ 

r?„(s,t)= V-r-Ti^^T^^— selK^tG [0,oo), (1) 

■ 1^ log{<^ {Zj{SnS,tnt))) 

converge weakly on the space of continous functions on M'^ x [0, oo) to the stationary Brown- 
Resnick process 

oo 

Vis,t) = \/ Cjexp{Wj{s,t)-6is,t)}, (2) 



where the deterministic function 6 is given in Assumption 2.1 and \ Wj{s, t),s £ M'^, t G [0, oo)}, 
j G N are independent replications of a Gaussian process with stationary increments, W{0, 0) = 
0, 'E{W{s,t)) = and covariance function for si, s 2 £W^,ti,t2 G [0,oo) 

Cov {W{si,ti), W{s2,t2)) = 6{si,ti) + 6{S2, t2) - 5{si - S2, ti - t2). 

The bivariate distribution function of rj can be expressed in closed form and is based on a well- 
known result by Hiisler and Reiss ]12^ : 

f 1 / log ^ \ 1 / loe; ^ \ 1 

F(xi,X2) = exp ^[ + --<!> Z^ + v^^(7^ , (3) 

where ^ denotes the distribution function of a standard normal distribution. 

Many correlation functions satisfy the following condition, which will be used throughout. 
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Condition 2.3. The correlation function has an expansion around zero, given by 



p{h, u) = 1 — Oi\\h\ 



, u G 



where < ai,a2 < 2 and 61,62 > 0. 



Condition 2.3 allows for an explicit expression of the limit function 6 in Assumption 2.1 

6{h,u) = 6i\\hr' +62\u 



|«2 



(4) 



where the scaling sequences (sn) and (tn) can be chosen as s„ — (logn)^/°i and t„ = (logn)^/"2_ 
The parameters 01,02 S (0,2] relate to the smoothness of the underlying Gaussian process in 
space and time, where the case oi = 02 = 2 corresponds to a mean-square differentiable process. 
For example, Gneiting's class of correlation functions [lOj satisfies Condition 2.3 For a detailed 
analysis of Gneiting's class and further examples we refer to Davis et al. [5j, Proposition 4.5, 
where the expansion around zero is calculated for several classes of correlation functions. A 
further property of the model defined in Proposition 2.2 is the closed form expression for the 
tail dependence coefficient, which is defined by 



where h = si — S2 denotes the spatial distance between two locations and n = ti — ^2 is the 
temporal lag. As derived in Section 3 in Davis et al. [5j, we obtain 



{h,u) = 2(1- <S>{^y5{h,u))) = 2(1- +02|'u|°2 



(5) 



3 Pairwise likelihood estimation 

In this section, we describe the pairwise likelihood estimation for the parameters of the model 
in ([1]) introduced in Section [2] Composite likelihood methods have been used, whenever the full 
likelihood is not available or intractable. We present the general definition of composite and 
pairwise likelihood functions in Section 3.1 Afterwards, we describe the details for our model. 



3.1 Basics on composite likelihood estimation 

Composite likelihood methods go back to Besag [1] and Lindsay [16] and there is vast literature 
available, from a theoretical and an applied point of view. For more information we refer to 
Varin p3] who gives an overview of existing models and inference including extensive references. 
In the most general setting, the composite log-likelihood function is given by 

lc(.il^,x) = ^Wilogf{x G Ai;'il)). 
1=1 

From this general form, special composite likelihood functions can be derived. For our setting 
we define the (weighted) pairwise log-likelihood function by 

n n 

PL{if);x) = ^^Wij\ogf^{xi,Xj), (6) 

2=1 j = l 
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where x = (xi, . . . , Xn) is the data vector, f,^{xi,Xj) is the density for the bivariate observations 
{xi,Xj) and Wij are weights which can be used for example to reduce the number of pairs 
included in the estimation. The parameter estimates are obtained by maximizing ([g]). 

As noted in Cox and Reid for dependent observations, estimates based on the composite 
likelihood need not be consistent or asymptotically normal. This is important for space-time 
applications, since all components may be highly dependent across space and time. 

3.2 Application to spatio-temporal max-stable random fields 

To derive pairwise-likelihood functions for the model defined in Proposition ([T]) we first need to 
derive the bivariate density function for the space-time max-stable process. For later purposes 
we state the closed form expression in the following lemma. Throughout we denote by $ and 
(/? the cumulative distribution function and the density of the standard normal distribution, 
respectively. For simplicity we suppress the argument {xi,X2)- 

Lemma 3.1. Set 6 := S{h, u) as given in Q and define for xi, X2 > 

V ^Mi'i') + (8) 

X\ X'2 

The partial derivatives of q'^ and q'^ are given by 

fl (1) o (1) o (2) a (2) 



dxi 2y/5xi 9x2 2-Jlx2 9x2 2^/5x2 dxi 2VSxi' 



The first and second order partial derivatives of V are given by 



dV 

dxi 
dV 
dx2 

dxidx2 



1 



X 



X. 



2^4^) 
1 



1 



1 



2V-Sxf^''i'^^ 2^1x^X2''^'^ 



(2)^ 



1 



2Vda 



1 



2yJ~5xi 



2^/5 



2V6 

45X1X2 



X2 

(2) 

^V- / {2)n 



Finally, the bivariate log-density is 

log /(xi, 3:2) = -V + log 



/'dV\ /dV\ _ d'^V 
\dxi J \dx2 J dxidx2 



(9) 



The resulting parameter vector is "0 = {9i,ai,92,ct2)- We first define the pairwise likelihood 
for a general setting with M locations si, . . . , sm and T time points < ti < • • • < < 00. In 
a second step we assume that the locations lie on a regular grid and that the time points are 
equidistant. 



M M T-l 

EE 

1=1 j=i+l k=l l=k+l 



EE 



wl I log f^{r]{si,tk),v{sj,ti)), 



(10) 
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where w\^^ > and w^^i' > denote spatial and temporal weights, respectively. Since it is 
expected that space-time pairs, which are far apart in space or in time, have only little influence 
on the dependence parameters to be estimated, we define the weights, such that in the estimation 
only pairs with a maximal spatio-temporal distance of {r,p) are included, i.e., 

'^IT = ^{\\s^-s,\\<r}^ = ^{\tu-ti\<p}, (11) 

where || • || denotes any arbitrary norm on R'^. The pairwise likelihood estimates are given by 
(^1,01,^2,02)= arg max PL(^^'^)(6'i, 01,^2,02). (12) 

(6*1, 01,62, 02) 



Using the definition of the weights in (11), the log-likelihood function in (10) can be rewritten 
as 

M M T-pmin{fc+p,T} 

PL(^^.^)(^) = ^ E E E ^ogU{i^{s.,,t,),i^{s,,ti)). 

1=1 j=i+i k=l l=k+l 
\\si-Sj\\<r 

The following sampling scheme is assumed throughout. 
Condition 3.2. We assume that the locations lie on a regular d- dimensional lattice, 

S = {sj,3 = 1, . . .,m'^] = = (n, . . .,id), ii, . . . ,Zd e {1, • • • }• 

Further assume that the time points are equidistant, 

= to < ti < ■ ■ ■ < tr < 00, \tk - tk-i\ = I, k = l,...,T. 



For later purposes, we rewrite the pairwise log-likelihood function under Condition 3.2 in the 
following way. Define Tir as the set of all vectors with non-negative integer-valued components 
h without the 0- vector, which point to other sites in the set of locations within distance r. Nott 
and Ryden [TTj call this the design mask. We denote by \7ir\ the cardinality of the set Tir- In 
our application, we will use the following design masks according to the Euclidean distance; 

7^1 = {(1,0), (0,1)} 

7^2 = ^iU {(1,1), (0,2), (2,0)} 

^3 = 1^2 U {(1, 2), (2, 1), (2, 2), (0, 3), (3, 0)} 

ni = n-iU {(1, 3), (3, 1), (2, 3), (3, 2), (4, 0), (0, 4)} 



Using Condition 3.2 and the design mask, the pairwise log-likelihood function in (12) can be 
rewritten as 



PL('"''^)(^) 



m m 



T min{A:+p,T} 

E'"E E E E ^''SU{v{sin,...,i,),tk),visih,...,i,)+h,ti)) 



«i=l *d=l henr k=l l=k+l 



m m T 



+hes 



E • • • E E ^V- in, .■■,id,k; nr,p) - 7^(- '^)(^), 

41=1 id=l k=l 



(13) 
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where 

k+p 

g^{ii,...,id,k;nr,p) = ^ ^ log /^(?7(s(i^^,.. i^), tfc), ry(s(j^_.,^j^) + h., t/)), (14) 

heHr l = k+l 

and 7^('^''^)('0) is a boundary term, given by 

mm, T k+p 

7^(-^^)(^) = E E E EE log/v,(r?(^(,„...,,),t,),7?(s(,,,...,,) + (15) 

S(^^,...,^^)+HS 1>T 




Figure 3.1: Visualization of the boundary term T?.*^™"^'"^^ for d = 2, m = 6 and any time point, where 
the set S of locations is the inner square and the outer polygon represents the endpoints of pairs in the 
boundary 



Figure 3.1 depcits a spatial grid with length m = 6, where the inner square is the set of 
observed locations S and the points in the outer polygon are endpoints of pairs which are in the 
boundary term TZ^"^ The figure visualizes the case which is represented by the quarter 
circles. 



4 Strong consistency of the pairwise likelihood estimates 

In this section we establish strong consistency for the pairwise likelihood estimates introduced 



in Section 3.2 For univariate time series models Davis and Yau [7j proved strong consistency of 
the composite likelihood estimates in full detail. For max-stable random fields with replicates, 
which are independent in time, Padoan et al. |T8] showed consistency and asymptotic normality 
for the pairwise likelihood estimates. In contrast to previous studies, where either the spatial 
or the time domain increases, we show strong consistency as the space-time domain increases 
jointly. 
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4.1 Ergodic properties for max-stable processes 

Stoev and Taqqu [21j introduced extremal integrals as an analogy to sum-stable integrals. Based 
on the extremal integral representation of max-stable processes Stoev |2D] establishes conditions 
under which the max-stable process is ergodic. Wang et al. [SH] extend these results to a spatial 
setting. In the following, let T(^hi,...,ha,u) denote the multiparameter shift-operator. In accordance 
with the definitions and results in Wang et al. [26], we define ergodic and mixing space-time 
processes. 

Definition 4.1. A strictly stationary space-time process [r]{s,t),s € 'K'^,t G [0,oo)} is called 
ergodic, if for all A,B e a{rj{s,t),s GR'^,t€ [0,oo)} 

mi m^i T 

hi=l h^=l u=l 

If the process satisfies additionally 

Jim P (a n ^,i,„,...,s,,„,t„){B)) = P{A)P{B), (17) 

for all sequences {(si^„, . . . , Sd,n, tn),n G N} with max{|si^„|, . . . , \sd,n\, \tn\} — ^ oo, we call the 
process mixing. 

Note in ( 16 ) that in contrast to the ergodic theorem in Wang et al. p6| , the number of terms 
in each sum is not equal, since we have an additional sum for the time component. We focus on 
max-stable processes with extremal integral representation 

e 

7]{si, ...,Sd,t)=J U^si,...,s,,t)if)dMi, (18) 

E 

where U(_si,...,sa,t) ■ ^Hm) ^H/^) given by U(_si,...,sa,t){f) = f ° T{si,...,sa,t) is a group of max- 
linear automorphisms with ?7(o,...,o,o)(/) = / a^^d the control measure fi is the distribution of the 
space-time process, 

/x(A) = P (7?(si, ...,Sd,t) e A), A£a |ry(s, t), {s, t) € M'^ x [0, cx))} . 

The following result is a direct extension of the uniparameter theorem established in Stoev |2; 
Theorem 3.4, and the multiparameter counterpart in Wang et al. [26|. 



Proposition 4.2 (Wang et al. |26j. Theorem 5.6). The max-stable process defined in (18) is 
mixing, if and only if 



t/(.i,„,...,.,,„,t„)(/) A C/(o,...AO)(/)rf^ = / f^(si,„,...,.,,„A)(/) A fdf^ ^ 0, (19) 

Je 

for all sequences {(si,„, . . ■,Sd,n,tn)} with max{|si,„|, . . . , \sd,n\, \tn\} oo as n oo. 

Wang et al. [26] showed, that the ergodic theorem stated above holds for mixing max-stable 



processes with extremal integral representation (18) in the case oiT = m. The extension to the 
multiparameter case where T 7^ m is a simple generalisation using Theorem 6.1.2 in Krengel 
[15j, which is a multiparameter extension of the Akcoglus ergodic theorem. Ergodic properties 
of Brown-Resnick processes have been studied for the uniparameter case in Stoev and Taqqu 
|21j and Wang and Stoev [27J. We summarize the results in the following proposition. 
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Proposition 4.3. The Brown- Resnick process in Proposition 2.2 with extremal integral repre- 
sentation 

e 

j Q:^v{W{s,t) - 5{s,t)]dMi s € M'^, t G [0, oo) 

E 

is mixing in space and time. The strong law of large numbers holds; 

^ rrt m T 

^ E • • • E "4- E [<7(r?(«(i,...,i),ti))] , mT ^ oo (20) 

ii=l id=l fc=l 

where g is a measurable function. 
4.2 Consistency for large mT 



In the following we show that the pairwise likelihood estimate resulting from maximizing (13) 



for the model defined in Proposition 2.2 is strongly consistent. 



Theorem 4.4. Assume that the true parameter vector ip* = (0J', a*, ^^^■^ compact 

set which does not contain and which satisfies for some c > 

^ C {min{0i,02} > c,ai,a2 G (0,2]}. (21) 
Assume also that the identifiability condition 

•0 = ^ <^ f,p{v{si,ti),ri{s2,t2)) = f^{ri{si,ti),ri{s2,t2)), a.s. (22) 
is satisfied for all (si,ti), {s2,t2). It then follows that the pairwise likelihood estimate 

ip^d^T = argmaxPL("^''^)(-0) 

is strongly consistent, i.e. '^^n'^ y ~^ V'* o,^ ~^ 00. 



Remark 4.5. For the identifiability assumption (22) we consider different cases according to 



the maximal space-time lag {r,p) included in the composite likelihood. Recall that the pairwise 



density, see Lemma 3.1, depends on the spatial distance h and the time lag u only through the 
function 5{h,u) = 9i\\h\\°'^ -\-62\u\"^. For specific combinations of {r,p) not all parameters are 
identifiable. For example, if the maximal spatial lag taken into account in the estimation equals 
one, i.e. r = 1, and p > 1, the parameter ai is not identifiable. Strong consistency still holds for 
the remaining parameters. Table [T] lists the various scenarios. 



Proof of Theorem 4-4 To show strong consistency of the estimates we follow the method of 



Wald P5|. Accordingly, it suffices to show the following conditions. 
(CI) Strong law of large numbers: Uniformly on the compact set ^, 

-^PL(-''^)(^) ^ PL(V') := E [5^(1, . . . , 1, 1; nr,p)] , mT ^ 00. 
(C2) The function PL{ip) is uniquely maximized at the true parameter vector -0* G ^. 
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Maximal spatial lag r 


Maximal temporal lag p 


Identifiable parameters 





1 


02 





p, p > 1 


O2, "2 


1 





6'i 


r, r > 1 





6*1, «! 


1 


1 


^1, ^2 


1 


p, p > 1 


01, 62, a2 


r, r > 1 


1 




r, r > 1 


p, p > 1 


Ol, «!, 6*2, 02 



Table 1: Identifiable parameters corresponding to different maximal space-time lags (r, p) included in 
the pairwise likelihood function. 



From (CI) and (C2) strong consistency follows. First we prove (CI). Recall from (13) that the 
pairwise likelihood function is given by 



m T 



ii=l «[i=l fc=l 



where and TZ^"^ (-0) are defined in ( 14 ) and ( 15 ) , respectively. The pointwise convergence of 



the first term on the right hand side to PL^ij^) follows immediately from Proposition 4.3 together 
with the fact that in (14) is a measurable function of lagged versions of rji^{ii,...,id),tk)- 



In the following, we show that the convergence is uniform and that the boundary term 



defined in ( 15 ) converges to zero almost surely. For both steps, observe first that we can bound 



the log-density from Lemma 3.1 For xi, X2 > 



|log/^(xi,X2) 



< 



-V + log 
1 



Xi 



dV dV 

dxi dx2 

+ 



dxidx2 



1 1 

< — + — + 



1 



Xl 



+ 



X2 
1 



+ 



1 



+ 



dV dV 



1-^2 



1 



dxi dx2 dxidx2 
1 1 1 



+ 



46{h,u) 



+ 



1-^2 



x\x2 



+ 



X1X2 



1-^2 



+ 



x\x2 



+ 



X?X^ 



+ 



1-^2 



X1X2 



+ 



xfx2 



+ 



2 



/y 2 ™2 

■^1-^2 



+ 



Jl) 



xfx2 



+ 



X1X2 



where q~^\ q^'^ and V are defined in ([T]) and ([s]), respectively, where <^(-) < 1 and ip{-) < 1 
were used. Since the marginal distributions of the max-stable space-time process are assumed 
to be standard Frechet, it follows that for every fixed location s G 5 and fixed time point t € T 
l/r]{s,t) is standard exponentially distributed. Using Holder's inequality, it follows that 



V [\^ogU{rj{si,ti),rj{s2,t2))\] <Ki + 



K2 



+ 



K3 



2^6{h,u) ' 4(5(^,n)' 

where Ki,K2, > are finite constants. Since the parameter space ^ is assumed to be compact 



and together with assumption (21), 5 can be bounded away from zero, i.e. 



5{h,u) = Oi\\hr^ +e2|ur' >min{0i,02}(||/i| 
>c(||h|r^ + |nr2) >c>0, 



+ \u 



02 > 



(23) 
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since 01,02 G (0,2], where c > is some constant independent of the parameters. Therefore, 



E, 



|log/v,(^(si,*i),^(s2,t2))|] <Ki + 



K2 ^ K3 

2Vd 4c 



= : < 00, 



(24) 



where K4 > 0. Note that in the same way we can show that the expectation of the squared 
bivariate log-density is finite, since it only involves higher order moments of the exponential 
distribution. 

To establish uniform convergence, we follow Straumann and Mikosch [22], Theorem 2.7, and 
show that 



E 



sup \g^{l,...,l,l;'Hr,p)\ 



< 00. 



It is sufficient to verify that 



E 



sup |log /v,(?/(s{n,...,id) , tk),visih,...,id) + h,tk + u)) \ 



< 00. 



Since the pairwise density is continous and because of the compact parameter space, the state- 



ment follows immediately using (23) and (24). 



As a last step for (CI) we show that the boundary term TZ^"^ converges to almost 

surely. For notational simplicity, we only consider the case d = 2. The general case is proved 



analogously. First note from (|15|) that 



E 



T k+p 



Yl IE [|log/^(r/(s(i^,i2),tfc),r/(s(i^,i^) +^,tO)| 



il = l 42 = 1 heHr fc = l i = fc+l 

S(,^.,2)+h^S 1>T 



< 



^ mm 
27^ X] X] X] 



i\=l 42=1 



S(r^,r^)+hiS 1>T 



yy K,< ^ 

^ ^ mT 

k=l l=k+l 



0, mT — >• cx). 



where we used the bound derived in (24) and the fact that the number of space-time points 



in the boundary is of order m (independent of T) and, therefore, can be bounded by K^m 
with > a constant independent of m and T. We write 7^^™ ''^H'lp) in the following way 



using the function in (14). Denote by Bm,T the set of space-time indices {ii,i2,k) for which 
+ h ^ S or k > T.ln 

Figure [3.1 1 i3m,T corresponds to the indices of the locations j^-j G S 
for which s^j^ j^) + ^ is in the outer polygon. The cardinality of the set Bm,T can be bounded by 
using the maximum norm instead of the euclidean norm, i.e. 

I'Sm.Tl < r(2m + 1) =: K^m. 

The boundary term in ( 15 ) can then be written as 

n^rn\T)^^^^ 5^(n,i2,A:;?^.,p), 



where is defined in (14). In the same way as before, it follows by the strong law of large 
numbers that 

1 



(«i,«2,fc)eB,„,T 



g^p{k,i2,k;'Hr,p) ^ E [g^{l, 1, l;nr,p)] , 
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uniformly on the compact set ^. Therefore 
1 



since E [|(7^(1, 1, 1; < oo. This proves (CI). 

To prove (C2), note by Jensen's inequahty that 



log 



/V;(X1,X2) 
/^*(X1,X2) 



< log ( E^. 



and, hence, 



PL{tp) < PL{ip*). 

So, "0* maximizes PL(ip) and is the unique optimum if and only if there is equality in Jensen's 



inequality. However, this is precluded by (22) 



□ 



5 Asymptotic normality of the pairwise likelihood estimates 

In order to prove asymptotic normality of the pairwise likelihood estimates resulting from max- 



imizing (13) we need the following results for the pairwise log-density. The proofs can be found 



in Appendix [A} 

Lemma 5.1. (1) The gradient of the bivariate log-density satisfies 



E^. 



IV^, \ogf^{r]{si,ti),r]{s2,t2))\' 



< oo 



(2) The Hessian of the pairwise log-density satisfies 



E, 



sup |V^log/^(r/(si,ti),r/(s2,t2) 



< OO. 



Assuming asymptotic normality of the pairwise score function, then it is relatively routine to 
show that the pairwise likelihood estimates are asymptotically normal. To formulate the result. 



recall from (13) that the pairwise likelihood function can be written as 

m m T 



il=l 



=1 fc=l 



where is defined in (14). The pairwise score function is then given by 

m m T 

E ■ ■ ■ E E ^V-^vCn, .■■,id,k; Ur,?) - Vv,7^('"'•^)('0), 

ii=l «(i=l k=l 

where V.^g^{ii, . . . ,i(i,k; Hr-.p) is the gradient of the function g^ with respect to tj). 

Theorem 5.2. Assume that the conditions of Theorem \4-4\ hold. In addition, assume that a 
central limit theorem holds for V^(7^(zi, . . . ,i(i,k; H-riP) in the following sense 

^ m m T 

= E •■■ E E^V'5v*(^i'---'^d'^;^'''f') — >AA(0,S), mT ^oo, (25) 



ii = l *d=l k=l 



12 



where ip* is the true parameter vector and S is some covariance matrix. Then it follows that the 
pairwise likelihood estimates Tp^nd x satisfy 

m'^/2^(-0„d^r - iP*) A AA(0, F-^^{F-Y), mT ^ oo, 

where 

F = E^. [-V2,5^.(l,...,l,l;-H.,p)] . 

Proof. We use a standard Taylor expansion of the pairwise score function around the true 
parameter vector: 

rW^/'^i^^d^r -r)=-( ^V^PL(™^^)(^))" (^^^V^PL(-^^)(^*)) 



m 



ii=l id=l k=l 



( m m T \ 

ii=l i2=l fc=l / 



where V £ ['V'm'* 'V'*]- For now, we ignore the boundary term and analyze the first terms in 
the outer brackets. By (25) the second term converges to a normal distribution with mean 
and covariance matrix S. For the first part we use the same arguments as in the consistency 
proof and show a strong law of large numbers. Since the underlying space-time process in the 
likelihood function is mixing, it follows that the process 

|v^(7^(ii, . . . , id, fe;'Hr;P), GZ'^,t/c G z| is mixing as a measurable function of mixing 

and lagged processes. To prove the uniform convergence we need to verify that 



sup |V^5^(l,...,l,l;^r,p)| 



< oo. 



This follows immediately from Lemma |5.1[ Putting this together with the fact that S 
[■j/j^d and because of the strong consistency, it follows that 



m T 



m 



^ E ■ ■ ■ E E^V'f^(^i' • • ■ ,id,k;nr,p) ^ E^* [V^5fv*(l> l;'^r,p)] =: -F. 

il = l id = ^ k=l 

For the boundary term TZ^"^'^''^^ observe that it can be written as 

{ii,...,ia,k)eBm,T 

Using assumption ( 25 ) together with the strong law of large numbers for | V^(7^ (ii , . . . , id, k; T-lr,p) | 
it follows in the same way as in the proof of Theorem 4.4 that 

Lv5,7^(-''^) ^ 0, and _l^Vv,7^(-''^) A 0. 



Combining these results, we obtain 



m 



m 



'^I^Vf{^^d T-tp*)^ AA(0, F'^J:{F~Y), mT ^ oo. 



□ 



In the next section we provide a sufficient condition for (25). 
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5.1 Asymptotic normality and «-mixing 

In this section we consider asymptotic normality of the parameters estimates for the process 
in Proposition |2.2[ Under the assumption of a-mixing on the random field the key is to show 
asymptotic normality for the score function of the pairwise likelihood. For an increasing time 
domain and fixed number of locations asymptotic normality of the pairwise likelihood estimates 
was shown in Huser and Davison [11]. The main difference between a temporal setting and a 
space-time setting is the definition of the a-mixing coefficients and the resulting assumptions to 
obtain a central limit theorem for the score function. 

We apply the central limit theorem for random fields established in Bolthausen [2] to the 
score function of the pairwise likelihood in our model. In a second step we verify the a-mixing 
conditions for the max-stable process introduced in Section [2j First, we define the a-mixing 
coefficients in a space-time setting as follows. Define the distances 

d{{si,ti), {s2, t2)) = max <^ max \si{i) - S2{i)\,\ti - ^2! L ■Si, ■§2 G Z'^, G N 

d(Ai, A2) = inf ti), (S2, ^2)), {siM) G Ai, (^2, ^2) e A2} , Ai, A2 C x N. 

Let further = a {r]{s,t), (s,t) € Aj} for i = 1,2. The mixing coefficients are defined for 
k,l,n > hy 

ak,iin) = sup n A2) - P{Ai)P{A2)\ : Ai G -Fa,, | Ai| < k, IA2I < I, d(Ai, A2) > n} (26) 

and depend on the sizes and the distance of the sets Ai and A2. A space-time process is 
called a-mixing, if afc ;(n) — )• as n — )• 00 for all k,l > 0. We assume that the process 
{r]{s,t),{s,t) G Z'^ X N} is a-mixing with mixing coefficients defined in (|26l), from which it 
follows that the score process 

{V^g^ih, ...,id,k; nr,p), tfc) G Z'^ x n} . (27) 

is a-mixing. We apply Bolthausen's central limit theorem to the process in (27). By adjust- 
ing the assumptions on the a-mixing coefficients to the score process, we obtain the following 
proposition. 

Proposition 5.3. Assume, that the following conditions hold: 

(1) The process \^{r]{s,t), {s,t) G Z^ x N| is strongly mixing with mixing coefficients ak^i{n) as 
in @. 

00 

(2) Z < 00 for k + l< 4(1^,1 + l){p + 1) and a(^\nr\+i){p+i),M = o{n-(^+^^). 

n=l 



(3) There exists /3 > such that 



E [|V^c/^.((ii, . . . ,id),k,nr,p)\ 



2+P 



< 00 and 



^^'^Oimr\ + l)ip+l)X\nr\ + l){p+l)in)^^^'^^^'^ < 

n=l 



Then, 



T 



^ • • • ^ ^ V^5^*(n, . . . ,id,k;Ur,Q) M{0, S), mT 00, 



where^= Yl J2 '^ov {V^g^-il, . . . ,l,l;nr,p),V^g.^*{ii, . . . ,id,k;nr,p)) . 
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In a second step, we want to analyze the strong mixing property and the related assumptions 
for our model. Recent work by Dombry and Eyi-Minko [t9] deals with strong mixing properties 
for max-stable random fields. By using point process representation of max-stable processes 
together with coupling techniques, they showed that the a-mixing coefficients can be bounded 
by a function of the tail dependence coefficient. A direct extension to the space-time setting 
gives the following lemma. 

Lemma 5.4 (Dombry and Eyi-Minko [9], Corollary 2.2). Consider the stationary max-stable 
space-time process {^^(s, t), (s, t) G Z*^ x N} with tail dependence coefficient x{h,u). The a- 
mixing coefficients in ( |26[ ) satisfy 

Oik,l{n) < kl sup x{h,u) and ak,oo{n') < k x{h,u). 

max{||^i||,|«|}>n ' max{||h||,|«|}>n 

For the model described in Proposition 2.2 with tail dependence coefficient x in ([s]), it follows 
by using the inequality for the normal tail probability ^{x) < e~^ that 



sup sup exp < - 

ll'^llil"!!^'^ max|||h|j,|M||>n L 



ak,i{n) < Akl 

max{||/i||,|'(t|}>n max{||h|j,|M|}; 

= m sup exp|-^(0i||^||°i + 02|u|°2; 



6{h, u) 



max{||h||,|M|}>n I ^ 

<m sup exp|-^min{^i,02}(max{||/i||,|n|})'"'"^"i'°2} 

max{||/i||,|ii|}>n L ^ 

For n — )• oo, this tends to zero for all k,l>0. Thus, {r]{s, t), (s, t) £ Z'^ x N} is strongly mixmg. 
This shows the first assertion in Corollary 5.3 Furthermore, for k + l < 4(|?^^| + 1) the 

coefficients satisfy 



^n'^ak,i{n) <m^n'^ sup exp i + 02|ur') i 

n=l ' n=l max{||;i||,|n|}>n I 2 J 

< 4A:/^n'^exp|-^min{ei,^2}ra'^'''^"i'"^^| < oo. 



n=l 

In addition, 

^'+'«(|W,.|+i)(p+i),oo(n) < n'^+^d^.l + l)5^exp|-^min{ei,02}x--^"^'-2> 



which proves (2) in Proposition 5.3 As for (3), from Lemma 5.1 and using /3 = 1 we know that 



E 



|Vv>5v*((^i> • • ■,id),k,nr,p) 



(2+/3) 



< OO. 



Using the same arguments as above, it is easy to see that the second condition in (3) holds. By 
combining the above results with Theorem 5.2 we obtain asymptotic normality for the parameter 
estimates i/'md t increasing number of space-time locations. 



Theorem 5.5. Assume that the conditions of Theorem 4-4 hold. Then, 

(m^T) _ 4 AA(0, F-iS(F-i)T), mT ^ oo. 
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with 
and 



S= ^ '^Cov{\/^g^*{l,...,l,l;'Hr,p),V^g^"{ii,...,id,k;nr,p)) . 



Si. 



6 Simulation study 

We illustrate the small sample behaviom' of the pairwise likelihood estimation for spatial dimen- 
sion d = 2 in a simulation experiment. The setup for this study is: 

1. The spatial locations consisted of a 10 x 10 grid 

^ = {^{h,i2) = {h,i2),k,i2 G {1, . . . , 10}} . 
The time points are chosen equidistantly, 1 < • • • T = 100. 

2. One hundred independent Gaussian space-time processes Zj{snS,tnt), j = 1, . . . , 100 were 
generated using the R-package RandomFields with covariance function p{snh,tnu). We 
use the following correlation function for the underlying Gaussian random field. 



p{h,u) = {1 + ei\\hr^ + e2\u 



a2\-3/2 



Assumption 2.1 is fuUfilled and the limit function 6 is given by 



lim logn(l - p{snh,tnu)) = 5{h,u) = ^e^Whr^ + ^02^1"'. 

3. The simulated processes were transformed to standard Prechet margins using the trans- 
formation -l/log(<I>(Zj(s,t))) for s e 5 and t E {t\, . . . .t-r) . 

4. The pointwise maximum of the transformed Gaussian random fields was computed and 
rescaled by 1/n to obtain an approximation of a max-stable random field, i.e. 

1 1 
' ' 100 ^.V log($(Z,(s„s,t„t)))' ' ^ ' 

5. The parameters 61,01,62 and Q2 for different combinations of maximal space-time lags 



{r,p) were estimated by maximizing (13). The program is adjusted such that it takes care 



of identifiability issues, when some of the parameters are not identifiable, cf. Remark 4.5 
6. Steps 1.-5. are repeated 100 times. 



Figures 6.1 and 6.2 show the resulting estimates as a function of {r,p), where the true 
parameter set is given by ■0* = {6^,01,62, 02) = (0.06, 1, 0.04, 1). Figure [6?T| shows the resulting 
estimates for the spatial parameters 61 and ai. The horizontal axis shows the different maximal 



space-time lags included in the pairwise likelihood function from ( 13 ). Each vertical dot shows the 
result for one specific simulation. The dotted lines show confidence bands based on the simulation 
results. In addition to the graphical output we calculate the root mean square error (RMSE) 
and the mean abolute error (MAE) to see how the choice of {r,p) influences the estimation. 
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We draw the following conclusions. As already pointed out by Davis and Yau |6j and Huser 
and Davison [11], there might be a loss in efficiency if too many pairs are included in the 
estimation. This can be explained by the fact that pairs get more and more independent as the 
space-time lag increases. Adding more and more pairs to the pairwise log-likelihood function 



can introduce some noise which descreases the efficiency. This is evident in Figure 6.2 for the 
temporal parameter a2, where the estimates vary more around the mean as more pairs are 
included in the estimation. 

An interesting observation for our model is that using a maximal spatial lag of or a maximal 
temporal lag 0, respectively, leads to very good results. For the spatial parameters, the space- 
time lags which lead to the lowest RMSE and MAE are (2, 0) for Oi and (2, 0) (RMSE) or (3, 0) 
(MAE) for 02 (see Table [2]), i.e. we use all pairs within a spatial distance of 2 or 3 at the same 
time point. Basically, this suggests that we could also estimate the spatial parameters based on 
each individual random field for fixed time points and then take the mean over all estimates in 
time. The same holds for the time parameters 62 and 02, where the best results in the sense of 
the lowest RMSE and MAE are obtained for the space-time lags (0, 3), i.e. if we use all pairwise 
densities corresponding to the space-time pairs (s, ti) and (s, ^2), where 1^2— ill < 3 (see Table[3|. 
The reason for this observation is that the parameters of the underlying space-time correlation 
function get "separated" in the extremal setting in the sense that for example a spatial lag equal 
to zero does not affect the temporal parameters Oi and ai and vice versa. 



ei 


(1,0) 


(1,1) 


(1,2) 


(1,3) 


(1,4) 


(1,5) 


(2,0) 


(2,1) 


(2,2) 


RMSE 


0.0123 


0.0118 


0.0121 


0.0122 


0.0123 


0.0124 


0.0103 


0.0104 


0.0105 


MAE 


0.0105 


0.0090 


0.0092 


0.0093 


0.0094 


0.0095 


0.0080 


0.0081 


0.0081 




(2,3) 


(2,4) 


(2,5) 


(3,0) 


(3,1) 


(3,2) 


(3,3) 


(3,4) 


(3,5) 


RMSE 


0.0104 


0.0104 


0.0104 


0.0106 


0.0107 


0.0108 


0.0108 


0.0107 


0.0108 


MAE 


0.0081 


0.0081 


0.0081 


0.0082 


0.0083 


0.0083 


0.0084 


0.0083 


0.0084 


ai 


(2,0) 


(2,1) 


(2,2) 


(2,3) 


(2,4) 


(2,5) 


(3,0) 


(3,1) 


(3,2) 


RMSE 


0.1338 


0.1398 


0.1530 


0.1492 


0.1543 


0.1569 


0.1351 


0.1409 


0.1579 


MAE 


0.1078 


0.1124 


0.1154 


0.1137 


0.1233 


0.1252 


0.1050 


0.1106 


0.1127 




(3,3) 


(3,4) 


(3,5) 


(4,0) 


(4,1) 


(4,2) 


(4,3) 


(4,4) 


(4,5) 


RMSE 


0.1596 


0.1639 


0.1649 


0.1423 


0.1483 


0.1614 


0.1673 


0.1735 


0.1751 


MAE 


0.1228 


0.1291 


0.1297 


0.1120 


0.1176 


0.1114 


0.1276 


0.1372 


0.1385 



Table 2: RMSE and MAE based on 100 simulations for the spatial estimates Ox and ai for different 
combinations of maximal space-time lags (r, p) . 
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(0,3) 


^n /i\ 
(0,4) 


(0,5) 


(1,2) 


^1 '^^ 
(1,3) 


(1,4) 


^1 K^ 
(1,5) 


(2,2) 


RMSE 


0.1317 


0.1269 


0.1280 


0.1289 


0.1442 


0.1401 


0.1426 


0.1438 


0.1463 


MAE 


0.1008 


0.0989 


0.1015 


0.1035 


0.1086 


0.1079 


0.1139 


0.1147 


0.1179 




(2,3) 


(2,4) 


(2,5) 


(3,2) 


(3,3) 


(3,4) 


(3,5) 


(4,2) 


(4,3) 


RMSE 


0.1532 


0.1580 


0.1619 


0.1473 


0.1531 


0.1589 


0.1642 


0.1549 


0.1607 


MAE 


0.1242 


0.1275 


0.1294 


0.1169 


0.1223 


0.1273 


0.1317 


0.1233 


0.1284 



Table 3: RMSE and MAE based on 100 simulations for the spatial estimates ^2 and 0:2 for different 
combinations of maximal space-time lags (r,p). 
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A Proof of Lemma 5.1 



In the following, we use the same abbreviations as in Lemma 3.1 The gradient of the bivariate 
log-density with respect to the parameter vector ip is given by 

^ . \ 51og/(xi,X2) 
Vv,log/(xi,X2) = 7^ V^6. 

Assume in the following that all parameters ^i, 01,^2 and 02 are identifiable. Since all partial 
derivatives 

a0\ 0O2 oai oa2 
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as well as all second order partial derivatives can be bounded from below and above for < 
min{||^||, |m|} ,max{||^||, |n|} < oo using assumption (21) and, independently of the parameters 
9i, 62, ai and 02, it suffices to show that 



E, 



51og/^(??(si,ti),ry(s2,t2)) 



and 



E, 



sup 



86 



92 logf^{r]{si,ti),r]{s2,t2)) 



36 



< 00 



< 00. 



Since 6 can be bounded away from zero using assumption (21), we can treat 5 as a constant. 
For simplification we drop the argument in the following equalities. Define 

dV ^ dV , , d^V 

Ai = t:—, A2 = T— , and A3 = . 

0x1 0x2 0x1X2 

The partial derivative of the bivariate log-density with respect to 6 has the following form 
dlogf^ _ dV 



85 



36 



+ {A1A2 - A3 



dA. . dA2 dA3\ 



36 ) 



We identify stepwise the "critical" terms, where "critical" means higher order terms of functions 
of xi and X2- To give an idea on how to handle the components in the derivatives, we describe 
one such step. Note that {A1A2 — ^43)"^ can be written as 



{A1A2 - As) 



-1 



X1X2 





1 


1 


1 




51 


X2 


X1X2 


^1 


i) 



where gi describes the sum of the components together with additional multiplicative factors. 
By using 

Mq^') and = ^^{qy] 



86 26 "^■'^ 
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where = log(x2/a;i)/(2\/5) + we have 



dAi 



Ao 



92 



(l)^2 



4^ (<0 



(l)^2 



•JU'^Jui2^ Jb-^Ju^J, Jb-^Jb'^ Jb-^Jb'^ Jb\^iJb<2 Jb^iJbi^ 



where §2 is a linear function of the components. By combining the two representations above, 

we obtain that all terms in 

{A1A2 - A'iY^{dAi/d6)A2 are of the form 



I logXil'^l] logX2|'^^ 
1 -^2 



, A;i,A;2,/c3, A:4 > 0. 



The second derivative of the bivariate log-density with respect to 6 is given by 



52 log d^v 



{d6f 



id6Y (^1^2-^3)- 



dA. dA2 dAs 

A2 + Ai 



d6 



86 86 



+ {A1A2 - A3 



( 8^Ai 
\{86r 



A2 + 2 



8A1 dA2 

1)6 1)6^ 



8^A2 8'^A3 
[86f {86f 



(28) 
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Stepwise calculation of the single components shows that all terms are also of form ( 28 ) . This 
implies that for both statements it suffices to show that for all /ci, ^2, ^3, ^4 > 



E 



(log7?(g,t))^i(logr?(s,t)) 
|r/(s,t)|'=3|r^(s,t)|fc4 



k2 



< OO. 



Since r]{s,t) is standard Frechet log(ry(s,t)) is standard Gumbel and l/r]{s,t) is standard expo- 
nential. Using Holder's inequality, we obtain 



E 



log(r?(g,t))|^i|log(r?(g,t))| 
\7j{s,t)\k3\r^(s,t)\''^ 



k2 



< E 



log{r]{s,t))\^^' E |log(r/(s,t))| 



4:k2 



1/2 



E 



E 



ri{s,t) 



4:k4 



1/2 



< OO, 



since all moments of the exponential and the Gumbel distributions are finite. 
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